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The Rayleigh equation ^R+RR+pp' 1 — with initial conditions -R(O) = Ro, R(0) = models the 
collapse of an empty spherical bubble of radius R(T) in an ideal, infinite liquid with far-field pressure 
p and density p. The solution for r = R/Ro as a function of time t = T/T c , where R{T C ) = 0, is 
independent of Ro, p, and p. While no closed- form expression for r(t) is known we find that 
ro(t) = (1 — t 2 ) 2 ^° approximates r(t) with an error below 1%. A systematic development in orders of 
t 2 further yields the 0.001%-approximation r*(t) = r (t)[l - a\ Li 2 . 2 i(i 2 )], where ai « -0.01832099 
is a constant and Li is the polylogarithm. The usefulness of these approximations is demonstrated 
by comparison to high-precision cavitation data obtained in microgravity. 
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I. INTRODUCTION 

George Gabriel Stokes might not have anticipated the 
importance of his endeavor when challenging his students 
in 1847 to calculate the collapse motion of an empty bub- 
ble in water [l| . The reach of this academic exercise was 
recognized in 1917 by Lord Rayleigh \2\ who conveyed a 
link between collapsing bubbles and the erosion damage 
found on ship propellers. The basic equation of motion 
of a collapsing bubble [3], now known as the Rayleigh 
equation (RE), reads 



dR 
dT 



d 2 R 
dT 2 



R+k = 0, 



(1) 



where the bubble radius R is a function of the time T, 
and k is a constant. Given the initial conditions 



R(0) = Ro 



dR 
dT 



= 



(2) 



T=0 



and the definition k = pp^ 1 , the RE describes the col- 
lapse of an empty spherical bubble of initial radius Rq 
in an incompressible, inviscid, infinite liquid with uni- 
form far-field pressure p and density p. If, alternatively, 
k is defined as (p — p Y )p^ 1 , the RE extends to the case 
of a gas- filled bubble with constant inner pressure p v . 
The RE neglects non-condens able bubble g well as 
surface tension and viscosity [Hj], liquid compressibility 
Q, and thermal effects @. However, regardless of those 
limitations and various enhanced models available today 
the RE remains widely used in practice, owing to 
its simplicity and often sufficient accuracy. This is de- 
spite the fact that the RE yields no closed-form solution 
for k > 0. While numerical solutions can be obtained, 
systematic analytical approximations offer better insight 
into the mathematical nature of the collapse, as we will 
show in this paper. Analytical approximations also be- 
come handy when the RE is integrated into multi-scale 
models, and they offer an intuitive understanding for the 
collapse motion. 



In this paper we develop highly accurate, yet remark- 
ably efficient analytical approximations for the solution 
of the RE. We first recall the standard normalization of 
the RE, based on which the analytical approximations are 
then obtained in a systematic way. These approximations 
are then compared against high-precision measurements 
of the most spherical bubbles available today. A short 
discussion concludes the paper. 



II. NORMALIZED RAYLEIGH MODEL 

Multypling Eq. (P) by 3R 2 (dR/dT), then integrating 
with respect to T, and expressing the integration constant 
using the initial conditions in Eq. @, we find 



dRY 
dTj 



^ ^ R" 



kR 3 = kR%. 



(3) 



Up to a factor Airp/3, Eq. Q expresses the conserva- 
tion of energy. Note that Eq. ([3]) only implies Eq. ([1} if 
R 2 (dR/dT) 7^ 0, and according to the initial conditions 
in Eq. ^ this relation breaks down at T = 0. Never- 
theless, we can still refer to Eq. ([3]) to analyze R(T) for 
T > 0. In particular, the collapse time T c , defined by 
R(T C ) = 0, is found by integrating dT from to T c and 
dR from Rq to 0. We obtain 

T c =£i?o£;- 1/2 , (4) 

where £ = ^/3/2 /^(r -3 - l)" 1 / 2 dr w 0.914681 is a uni- 
versal constant called the Rayleigh factor. Normalizing 
the radius to r = R/Rq and the time to t = T/T c , Eqs. ([T]) 
and ([2]) become 



2 
r(0) 



r 2 + rr + i 2 = 0, 



1, r(0) = 0. 



where dots denote derivatives with respect to t. 
the same normalization, Eq. @ translates to 



(6) 
Using 



(7) 
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Substituting Eq. ([7]) back into Eq. ([5]) then implies 

r = -er-\ (8) 

which is an interesting, but little known alternative form 
of the normalized RE given in Eq. (JSJ) - 

The key advantage of Eq. ([5]) over Eq. ([TJ is its invari- 
ance with respect to Rq and k. Stated differently, we only 
need to solve Eq. ([5]) once in order to solve Eq. ((T|) for 
any choice of Ro > and k > 0. The solution of Eq. (0 
in the range t = [0, 1] is displayed in Fig. QJb). This so- 
lution was obtained using a Cash-Karp fourth-fifth order 
Runge-Kutta method [1CJ- The relative error made on 
the collapse time lies below 10 -15 , and thus our numeri- 
cal solution for r(t) can be considered exact as far as this 
article is concerned. 



III. ANALYTICAL APPROXIMATIONS 

To hnd analytical approximations for r(t) we first no- 
tice that the Rayleigh model given in Eqs. (J5J) and ([5]) 
is symmetric in time. It follows that r(t) — r(— t) can 
necessarily be expressed as a function of t 2 . Second, we 
observe that r(t) is non-analytic at t = ±1, as can be 
seen from the divergence of r as r — > in Eq. ([7])- Since 
r(t) converges as t —> ±1, the singularities are neither 
essential singularities nor poles, and therefore must be 
branch points induced by the double-valued square root 
appearing when Eq. ([7} is solved for r. The simplest func- 
tion exhibiting such a pair of branch points at t = ±1 is 
the power law ro(t) = (1 — t 2 ) a with a <E ]0, 1[. In order 
to use ro(t) as an approximation of r(t) the parameter a 
can be determined in two ways. First, we can request that 
p'o(0) = r(0), which together with Eqs. © and ([5]) im- 
plies a = £ 2 / 2 ~ 0.418321. Alternatively, we can impose 
that ro(t) and r(t) exhibit similar asymptotic behavior at 

t = 1. In fact, Eq. §§§ implies r cx r^ 1 ^ and Eq. ([TJ 
implies r oc r~ 3 / 2 as t — > 1. Matching the powers of the 
two asymptotic functions we find a = 2/5 = 0.4. The 
coincidental close similarity of a determined at t = and 
t = 1 suggests that ro(t) is a good approximation of r(t) 
for both values of a. Here we choose a = 2/5 and the 
corresponding approximation 



ro(t)^(l-t 2 Y 



(9) 



hence accepting that r'o(0) ^ f(0) for the moment. This 
approximation is displayed in Fig. [ljb) . Its similarity to 
the Rayleigh function r(t) is demonstrated by the small 
residual ro(t) — r(t), shown in Fig. [UJc). Indeed, r*o(i) 
never differs by more than 0.01 from r(t), and by con- 
struction we have ro(t) = r(t) at t € {—1,0,1}. The 
normalized velocity residual [ro(i) — r(t)]/f(t), shown in 
Fig.[TJd), takes absolute values up to 0.044. The fact that 
limt-K) fa (i) — r(t)]/r(t) = f'o(0)/r(0) — 1 differs from zero 
explicitly reveals that f (0) =^ ^(0). 

We now improve the accuracy of our approximation 
ro(i) at t = through the modified ansatz ?'oo(i) = 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 
Normalized time t 

FIG. 1: (Color online) (a) Four subsequent high-speed images 
of the collapsing spherical bubble in microgravity. (b) Col- 
lapse functions: exact solution r(t) of Eqs. ([5]) and (JSJl (solid 
line), measurement r ^ B (t) (dots), and zeroth-order analytical 
approximation ro(t) (dashed line), (c) Errors of r b s (t) and 
the analytical models r n (t) [Eq. [11] and r* (t) [Eq. [13] relative 
to the exact solution r(t). Bars represent 67% statistical mea- 
surement uncertainties, (d) Errors of the velocities r b s (t), 
f,(t), and r n (t) relative to r(t). 
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TABLE I: Coefficients a q in the summations of Eqs. (|10[) and 
and accuracies of the approximations r q (t). 
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FIG. 2: Dots: analytically calculated coefficients a q in the 
expansion of Eq. (| 10p . Solid line: power law fit of Eq. (|12l) . 
which has been forced to pass through aj.. 



J"o (£)/(£), where /(£) is a smooth function defined by the 
condition that d 9 r oo (0)/dt 9 = d q r(Q)/dt q for all deriva- 
tives of order q > 1. As we will see, this infinite set 
of boundary conditions can be met by restricting /(£) 
to functions that can be expressed as Taylor series on the 
closed interval £ € [—1,1]. To respect the time-symmetry, 
/(£) must be even, i.e., all odd powers in the Taylor series 
vanish. Thus, 



,(£)^(l-£ 2 



9 =0 



/ j a qt 2<1 i 



(10) 



where a q are real constants 



The condition roo(0) = 
r(0) = 1 immediately implies clq = 1. All other coeffi- 
cients a q are obtained by matching the even-order deriva- 
tives of roo(£) and r(£) at t = 0. Those derivatives 
are evaluated analytically by differentiating Eqs. ([TO]) 
and ©, respectively, and applying the initial condi- 
tions given in Eq. ©. We can proceed iteratively: first 
use d 2<? r oo (0)/d£ 29 = d 2 «r(0)/d£ 2<? with q = 1 to get 
oi = 2/5 - £ 2 /2 ~ -0.01832099, then with q = 2 to 
get a 2 = 3/25 + 2oi/5 - £ 4 /6 « -0.00399003, and so 
forth. Note that all odd-order derivatives, such as r(0) 
and ^00(0), vanish due to time-symmetry. The analytical 
expressions for even-order derivatives of r(£) and ?'oc(£) 
get cumbersome as q increases; however, the coefficients 
a q are easily obtained by using analytical software tools. 
The numerical values of a q up to q = 10 are given in Ta- 



ble U and those up to q 
these values of t 
approximations 



22 are plotted in Fig. [2] Using 



these values of the coefficients a q we can now construct 



29 



(11) 



of any order n. Since a q < and \a q \ < |a.g-i| for all 
q > 0, the approximations r n (£) converge monotonically 
towards roo(£) as n — > 00. Figures QJc) and QJd) show 
that r n (t) and r n (t) numerically converge towards the 



Rayleigh solution r(£) and r(£). Already the second-order 
approximation r 2 (£) is roughly 10 times better than the 
zeroth-order approximation ro(£) discussed before. Yet, 
from a mathematical point of view, the crucial question is 
whether or not roo(£) is identical to r(£) for all £ G [—1,1], 
i.e. if r 00 (£) is the solution of the normalized RE. The 
answer is no, since r(t) in Eq. © and r oc (£) derived from 
Eq. (TTOj) do not obey the same asymptotic behavior as 
£ — > 1. In other words, r^it) remains an approximation 
of r(£), no matter the choice of the real coefficients a q . 

Figure [2] uncovers that q and a q exhibit a remarkably 
tight power-law relation. Although this relation is not 
analytically exact, we can approximate the values of a q 
for q > as 



a-i q 



-2.21 



(12) 



with the exact a x = 2/5 - £ 2 /2 « -0.01832099 given in 
Table H] The relation given in Eq. (|12p is plotted as a solid 
line in Fig. [2j Upon approximating the coefficients a q in 
7*00 (£) for q > by Eq. (fT2"j) we obtain 



(£)= (l-£ 2 ) 5 l + aiLi 2 . 21 (£ 2 ) 



(13) 



where Li s (x) = X)«*Li q~ s x 9 is the polylogarithm, also 
known as Jonquiere's function. Many programming lan- 
guages contain Li s (a;) in their standard libraries. We em- 
phasize that r*(£) slightly differs from r^t) since the 
latter uses the exact coefficients a q rather than those ap- 
proximated by Eq. (|12p . Nonetheless, r*(£) is a very pre- 
cise approximation of the Rayleigh solution r(£) as can 
be seen from the residuals (multiplied by a factor 100) 
shown in Figs. [TJc) andQJd). 

To quantify the accuracy of our approximations in a 
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more refined way we consider the measures 



e(r x ) = max \r x (t) - r(t) \ , (14) 




Here, e(r x ) is the maximal error and cr(r x ) the standard 
deviation of the approximation r x (t) with respect to the 
exact solution r(t), and likewise e(r x ) and &(r x ) give the 
accuracy of the derivative f x (t) . The residuals r x (t) — r{t) 
in Eqs. (TToT) and (JTTJ) are normalized relative to f(t) to 
obtain meaningful, converging measures as |f| — ¥ oo. It 
turns out that for all our approximations except tq (t) the 
absolute value of [r x {t)— f{t)]/r(t) is maximal as t — > 1. In 
this limit, the numerical evaluation of e is delicate because 
of the divergence of r x (t) and r(t). However, we can use 
the analytical expression lim t ^i = [r x (t) — f(t)]/f(t) = 
{V6a/£) a /(l) - 1, where f(t) = r x (t)/(l - t 2 ) 2 ? 5 is the 
sum on the right-hand sides of Eqs. (fTU)) , (fl~T1) . and (fT3")) , 
respectively. 

The logarithms of e, a, e, and a for various approxima- 
tions are given in Table [I] These values support and ex- 
tend the previous discussion of the residuals in Figs. [He) 
and [TJd). In particular, we find that r»(i) approxi- 
mates r(t) with an error below 10~ 5 on the whole interval 
t € [0, 1], while r*(i) approximates f{t) at a relative error 
below than 10~ 4 . 



IV. COMPARISON TO OBSERVED DATA 

We now present a state-of-the-art experiment (see de- 
tails in Ref. |9| ) of millimeter-sized bubbles with almost 
perfect sphericity. The hig h validity of the Rayleigh 
model for these bubbles [1J| stresses the need for accu- 
rate approximations such as r*(i) when analyzing such 
bubbles. 

In the experiment, single bubbles grow inside a liq- 
uid from a point=-plasma generated by a mirror-focused 
nanosecond laser pulse. The bubbles are sufficiently 
spherical that the hydrostatic pressure gradient caused 
by gravity becomes the dominant source of asymmetry 
in the collapse and rebound of the bubbles (see Fig. 1 (a) 
in Ref. @). To avoid this source of asymmetry the ex- 
periment is performed in micro-gravity conditions (ESA, 
53rd parabolic flight campaign). Therefore, the experi- 
ment can be considered as producing the most spherical 
cavitation bubbles available at present. 



The spherical bubble considered here has a maximal 
radius Rq = (2.786 ± 0.007) mm and collapses within 
T c = (508.28±0.10) [is under the driving pressure p—p v = 
(25.34 ± 0.15) kPa. This bubble is centered inside a vol- 
ume (178 x 178 x 150) mm 3 of demineralized water at 
(26 ± 0.5)°C. The bubble radius R ohs {T) is measured at 
sub-micron precision from a movie obtained with a high- 
speed camera, operating at inter-frame spacings of 10 /xs 
with exposure times of 370 ns. The high-s pee d movie and 
complementary data are available online [13j. 

Four selected time-frames of the collapsing bubble are 
shown in Fig. HJa). The evolution of the observed nor- 
malized radius r b s = Robs/Ro is plotted in Fig. Hfb). 
We find that r b s (t) closely follows the Rayleigh solution 
r(t), as emphasized by the residual r bs(i) — r(i), shown 
in Fig. []Jc). At no time does r ^ s (t) differ by more than 
10~ 3 from r(t). Therefore, considering the different ap- 
proximations r*(t) and r n (t) with n < 5, we see that only 
the accuracy of r*(i) is sufficient to compare the experi- 
mental data against the Rayleigh model. Such a compar- 
ison makes it possible, for instance, to efficiently analyze 
the remaining oscillatory residual r b s (t) — r*(i), which 
may be explained in terms of surface tension and vis- 
cosity [l3|, liquid compressibility @ and thermal effects 
[6]. Basic estimates of these effects unveil that the ex- 
cellent match between r Q b s (i) and r(t) is partially due to 
a compensation of surface tension, which accelerates the 
collapse, by compressibility and non-condensable gases. 



V. CONCLUSIONS 

In summary, we have developed analytical approx- 
imations to the solution r(t) of the RE, which pre- 
serve the time-symmetry, the boundary conditions at 
t = (to arbitrary order), the branch point singulari- 
ties at t = ±1, and the asymptotic behavior of r(t) at 
t = ±1. Despite their elementary forms, the zeroth- 
order approximation ro(t) yields a maximal error e(ro) of 
only about 0.01, whereas the best approximation r„(t), 
expressed in terms of the polylogarithm, reduces this 
error to below 10 -5 . Moreover, approximations r n (t) 
with any smaller accuracy can be systematically con- 
structed. For example, if 0.1% accuracy is desired then 
r 2 (t) « (1 - t 2 ) 2 / 5 (l - 0.01832£ 2 - 0.00399* 4 ) suffices. A 
comparison of these approximations against state-of-the- 
art measurements of highly spherical cavitation bubbles 
revealed that the residuals r»(t) — r(t) are more than 100 
times smaller than the observed residuals r ^ s (t) — r(t), 
shown in Fig. [He). Thus the approximation r*(t) is by 
far sufficient for all practical purposes. 

This work was supported by the Swiss NSF (200020- 
116641). 
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